function [e_Y_cpi,m_e_mj_cpi,cov_e_mj_cpi,check_mj_cpi,...
    e_Y_doubledef,m_e_mj_doubledef,cov_e_mj_doubledef,check_mj_doubledef,...
    e_Y_def,m_e_mj_def,cov_e_mj_def,check_mj_def] ...
    = elasticity_decomp_sectorCES_fun(delta,VA_hat_mj,...
    P_hat_n,RGDP_cpi_hat_n,VA_mj_0,GDP_deflator_n,RGDP_def_hat_n,DoubleDeflation_deflator_n,RGDP_doubledef_hat_n)
        

global shockT alphaT rhoT lambdaT etaT rho sigma lambda eta varphi ...
    tol maxit D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    N J FI_J france countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO sector_algorithm ... 
    sL_n sPi_n sD_n alpha_WIOD_j alpha_WIOD_francej labour_share_PWT

% Calculates the elasticity decomposition at the sector (ij) level

% Decomposition of aggregate real GDP elasticity into sector mean and
% covariance components

% Reshape VA_hat_mj so that it is country*sector and then deflate by
% country change in inflation to get real measures

%% Calculation using CPI

VV_mj = reshape(VA_hat_mj,[J,N])';
real_VA_hat_mj = VV_mj./repmat(P_hat_n',1,J);
clear VV_mj;

e_mj_cpi = (real_VA_hat_mj-1)/delta;
e_Y_cpi = (RGDP_cpi_hat_n-1)/delta;

m_e_mj_cpi = nanmean(e_mj_cpi,2) ; 
VV_mj = reshape(VA_mj_0, [J N])';
w_mj = VV_mj./repmat(sum(VV_mj,2),1,J);
w_mj_bar = size(w_mj,2)-1;% N - 1 correction for small small as in COV
clear  VV_mj;


cov_e_mj_cpi = zeros(N,1);
for m = 1:N
    C = nancov(w_mj(m,:)*w_mj_bar,e_mj_cpi(m,:));
    cov_e_mj_cpi(m) = C(1,2);
    clear C;
end

check_mj_cpi = e_Y_cpi-m_e_mj_cpi-cov_e_mj_cpi;

%% Calculation using GDP Deflator

VV_mj = reshape(VA_hat_mj,[J,N])';
real_VA_hat_mj_def = VV_mj./repmat(GDP_deflator_n,1,J);
clear VV_mj;

e_mj_def = (real_VA_hat_mj_def-1)/delta;
e_Y_def = (RGDP_def_hat_n-1)/delta;

m_e_mj_def = nanmean(e_mj_def,2) ;
VV_mj = reshape(VA_mj_0, [J N])';
w_mj = VV_mj./repmat(sum(VV_mj,2),1,J);
w_mj_bar = size(w_mj,2)- 1; % N - 1 correction for small small as in COV

clear  VV_mj;

cov_e_mj_def = zeros(N,1);
for m = 1:N
    C_deflator = nancov(w_mj(m,:)*w_mj_bar,e_mj_def(m,:));
    cov_e_mj_def(m) = C_deflator(1,2);
    clear C_deflator;
end

check_mj_def = e_Y_def-m_e_mj_def-cov_e_mj_def;

%% Calculation using double-GDP Deflator

VV_mj = reshape(VA_hat_mj,[J,N])';
real_VA_hat_mj_doubledef = VV_mj./repmat(DoubleDeflation_deflator_n,1,J);
clear VV_mj;

e_mj_doubledef = (real_VA_hat_mj_doubledef-1)/delta;
e_Y_doubledef = (RGDP_doubledef_hat_n-1)/delta;

m_e_mj_doubledef = nanmean(e_mj_doubledef,2) ;
VV_mj = reshape(VA_mj_0, [J N])';
w_mj = VV_mj./repmat(sum(VV_mj,2),1,J);
w_mj_bar = size(w_mj,2)- 1; % N - 1 correction for small small as in COV

clear  VV_mj;

cov_e_mj_doubledef = zeros(N,1);
for m = 1:N
    C_doubledeflator = nancov(w_mj(m,:)*w_mj_bar,e_mj_doubledef(m,:));
    cov_e_mj_doubledef(m) = C_doubledeflator(1,2);
    clear C_doubledeflator;
end

check_mj_doubledef = e_Y_doubledef-m_e_mj_doubledef-cov_e_mj_doubledef;


end

